In this part of the analysis we apply Revelio algorithm to explore cell cycle dynamic of pallial and hem domain radial glial cells
library(Seurat)
library(Revelio)
library(princurve)
library(orthologsBioMART)
library(Matrix)
library(dplyr)
library(RColorBrewer)
library(ggplot2)
library(ggExtra)
library(cowplot)
library(wesanderson)
#Set ggplot theme as classic
theme_set(theme_classic())Hem.data <- readRDS("../QC.filtered.cells.RDS")DimPlot(object = Hem.data,
group.by = "Cell_ident",
reduction = "spring",
cols = c("#83c3b8", #"ChP"
"#009fda", #"ChP_progenitors"
"#68b041", #"Dorso-Medial_pallium"
"#e46b6b", #"Hem"
"#e3c148", #"Medial_pallium"
"#b7d174", #2
"grey40", #4
"black", #5
"#3e69ac" #"Thalamic_eminence"
))Idents(Hem.data) <- Hem.data$Cell_identChP.data <- subset(Hem.data, idents = c("ChP", "ChP_progenitors"))
DimPlot(ChP.data,
reduction = "spring",
pt.size = 1,
cols = c("#83c3b8", "#009fda")) + NoAxes()FeaturePlot(object = ChP.data ,
features = c("Fgf8", "Fgf17", "Adamts15", "Fgfbp1"),
pt.size = 0.5,
cols = c("grey90", brewer.pal(9,"YlGnBu")),
reduction = "spring",
order = T) & NoAxes() & NoLegend()ChP.data <- AddModuleScore(ChP.data,
features = list(c("Fgf8", "Fgf17", "Adamts15", "Fgfbp1")),
ctrl = 10,
name = "Septum")
FeaturePlot(object = ChP.data ,
features = c("Septum1"),
pt.size = 0.5,
cols = rev(brewer.pal(10,"Spectral")),
reduction = "spring",
order = T) & NoAxes()hist(ChP.data$Septum1, breaks = 20)ChP.data$Septal.prog <- ChP.data$Septum1 > 0.1p1 <- DimPlot(ChP.data,
reduction = "spring",
group.by = "Septal.prog",
pt.size = 1) + NoAxes()
p2 <- FeaturePlot(object = ChP.data ,
features = c("Fgf17"),
pt.size = 0.5,
cols = c("grey90", brewer.pal(9,"YlGnBu")),
reduction = "spring",
order = T) & NoAxes() & NoLegend()
p1 + p2ChP.data <- subset(ChP.data,
subset = Septal.prog == FALSE & ChP.data$Spring_1 > 1300)DimPlot(ChP.data,
reduction = "spring",
pt.size = 1,
cols = c("#83c3b8", "#009fda")) + NoAxes() ## Fit principal curve
Trajectories.ChP <- ChP.data@meta.data %>%
select("Barcodes", "nUMI", "Spring_1", "Spring_2")fit <- principal_curve(as.matrix(Trajectories.ChP[,c("Spring_1", "Spring_2")]),
smoother='lowess',
trace=TRUE,
f = 0.8,
stretch=2)## Starting curve---distance^2: 137968435096
## Iteration 1---distance^2: 49449688
## Iteration 2---distance^2: 45168261
## Iteration 3---distance^2: 44040505
## Iteration 4---distance^2: 43717320
## Iteration 5---distance^2: 43609264
## Iteration 6---distance^2: 43575985
#The principal curve smoothed
ChP.pc.line <- as.data.frame(fit$s[order(fit$lambda),])
#Pseudotime score
Trajectories.ChP$Pseudotime <- fit$lambda/max(fit$lambda)
#Inverse the score if positive correlation with progenitor marker
if (cor(Trajectories.ChP$Pseudotime, ChP.data@assays$SCT@data['Hmga2', Trajectories.ChP$Barcodes]) > 0) {
Trajectories.ChP$Pseudotime <- -(Trajectories.ChP$Pseudotime - max(Trajectories.ChP$Pseudotime))
}
ChP.data$Pseudotime <- Trajectories.ChP$PseudotimeFeaturePlot(object = ChP.data,
features = "Pseudotime",
pt.size = 2,
cols = rev(colorRampPalette(brewer.pal(n =10, name = "Spectral"))(100)),
reduction = "spring",
order = T) & NoAxes()Prog.data <- subset(ChP.data, idents = c("ChP_progenitors"))
DimPlot(Prog.data,
reduction = "spring",
pt.size = 1,
cols = c("#009fda")) + NoAxes()Prog.data <- NormalizeData(Prog.data, normalization.method = "LogNormalize", scale.factor = 10000, assay = "RNA")Cellcyclegenes <- revelioTestData_cyclicGenes
head(Cellcyclegenes)## # A tibble: 6 × 5
## G1.S S G2 G2.M M.G1
## <fct> <fct> <fct> <fct> <fct>
## 1 ABCA7 ABCC2 ALKBH1 ADH4 AFAP1
## 2 ACD ABCC5 ANLN AHI1 AGFG1
## 3 ACYP1 ABHD10 AP3D1 AKIRIN2 AGPAT3
## 4 ADAMTS1 ACPP ARHGAP11B ANKRD40 AKAP13
## 5 ADCK2 ADAM22 ARHGAP19 ANLN AMD1
## 6 ADCY6 ANKRD18A ARL4A ANP32B ANP32E
We use orthologsBioMART library to map human to mouse mouse orthologs
G1.S <- findOrthologsHsMm(from_filters = "hgnc_symbol",
from_values = as.character(Cellcyclegenes$G1.S),
to_attributes = "external_gene_name")
S <- findOrthologsHsMm(from_filters = "hgnc_symbol",
from_values = as.character(Cellcyclegenes$S),
to_attributes = "external_gene_name")
G2 <- findOrthologsHsMm(from_filters = "hgnc_symbol",
from_values = as.character(Cellcyclegenes$G2),
to_attributes = "external_gene_name")
G2.M <- findOrthologsHsMm(from_filters = "hgnc_symbol",
from_values = as.character(Cellcyclegenes$G2.M),
to_attributes = "external_gene_name")
M.G1 <- findOrthologsHsMm(from_filters = "hgnc_symbol",
from_values = as.character(Cellcyclegenes$M.G1),
to_attributes = "external_gene_name")
gene.list <- list(G1.S$external_gene_name,
S$external_gene_name,
G2$external_gene_name,
G2.M$external_gene_name,
M.G1$external_gene_name)
CCgenes <- t(plyr::ldply(gene.list, rbind))
colnames(CCgenes) <- colnames(Cellcyclegenes)rawCounts <- as.matrix(Prog.data[["RNA"]]@counts)# Filter genes expressed by less than 10 cells
num.cells <- Matrix::rowSums(rawCounts > 0)
genes.use <- names(x = num.cells[which(x = num.cells >= 10)])
rawCounts <- rawCounts[genes.use, ]rm(list = ls()[!ls() %in% c("rawCounts", "CCgenes", "ChP.data", "Prog.data")])
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 5781851 308.8 10748997 574.1 8259915 441.2
## Vcells 119685790 913.2 611435808 4664.9 621902469 4744.8
We can now follow the tutorial form the package github page
myData <- createRevelioObject(rawData = rawCounts,
cyclicGenes = CCgenes,
lowernGeneCutoff = 0,
uppernUMICutoff = Inf,
ccPhaseAssignBasedOnIndividualBatches = F)## 2021-12-09 17:27:00: reading data: 2.32secs
rm("rawCounts")
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 5783930 308.9 10748997 574.1 8259915 441.2
## Vcells 119738779 913.6 489148647 3732.0 621902469 4744.8
The getCellCyclePhaseAssignInformation filter “outlier” cells for cell cycle phase assignation. We modified the function to keep all cells as we observed this does not affect the global cell cycle fitting procedure
source("../Functions/functions_InitializationCCPhaseAssignFiltering.R")
myData <- getCellCyclePhaseAssign_allcells(myData)## 2021-12-09 17:27:07: assigning cell cycle phases: 13.65secs
myData <- getPCAData(dataList = myData)## 2021-12-09 17:27:29: calculating PCA: 10.25secs
myData <- getOptimalRotation(dataList = myData)## 2021-12-09 17:27:56: calculating optimal rotation: 1.47secs
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 5813634 310.5 10748997 574.1 10748997 574.1
## Vcells 161090841 1229.1 489148647 3732.0 621902469 4744.8
To obtain the best linear ordering of cell along G1 to M phase we compare the output `ccPercentageUniformlySpaced’ form Revelio and the one obtained with the Angle method describe in Saelens et. al. Nature Biotechnology 2019
New dataframe
CellCycledata <- cbind(as.data.frame(t(myData@transformedData$dc$data[1:2,])),
nUMI= myData@cellInfo$nUMI,
Revelio.phase = factor(myData@cellInfo$ccPhase, levels = c("G1.S", "S", "G2", "G2.M", "M.G1")),
Revelio.cc= myData@cellInfo$ccPercentageUniformlySpaced,
Seurat.cc= Prog.data@meta.data[myData@cellInfo$cellID,"CC.Difference"])Compute cycling trajectory with Angle
CellCycledata$Angle.cc <- atan2(CellCycledata$DC1, CellCycledata$DC2) / 2 / pi + .5ggplot(CellCycledata, aes(DC1, DC2)) +
geom_point(aes(color = Revelio.phase)) +
scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5]))p1 <- ggplot(CellCycledata, aes(DC1, DC2)) +
geom_point(aes(color = Revelio.phase)) +
scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5]))
p2 <- ggplot(CellCycledata, aes(DC1, DC2)) +
geom_point(aes(color=Seurat.cc), size=2, shape=16) +
scale_color_gradientn(colours=rev(colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)),
name='Seurat_cc')
p3 <- ggplot(CellCycledata, aes(DC1, DC2)) +
geom_point(aes(color=Angle.cc), size=2, shape=16) +
scale_color_gradientn(colours=rev(colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)),
name='Angle.cc')
p4 <- ggplot(CellCycledata, aes(DC1, DC2)) +
geom_point(aes(color=Revelio.cc), size=2, shape=16) +
scale_color_gradientn(colours=rev(colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)),
name='Revelio_cc')
(p1+p2)/(p3+p4)ggplot(CellCycledata, aes(x= Revelio.cc, y= nUMI/10000)) +
geom_point(aes(color= Revelio.phase), size=0.5) +
scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) +
geom_smooth(method="loess", n= 50, fill="grey") +
ylim(0,NA)Prog.data$Revelio.DC1 <- CellCycledata$DC1
Prog.data$Revelio.DC2 <- CellCycledata$DC2
Prog.data$Revelio.phase <- CellCycledata$Revelio.phase
Prog.data$Angle.cc <- CellCycledata$Angle.cc
Prog.data$Revelio.cc <- CellCycledata$Revelio.ccp1 <- FeaturePlot(object = Prog.data,
features = "Revelio.cc",
pt.size = 1,
cols = rev(brewer.pal(10,"Spectral")),
reduction = "spring",
order = T) & NoAxes()
p2 <- DimPlot(object = Prog.data,
group.by = "Revelio.phase",
pt.size = 1,
reduction = "spring",
cols = c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) & NoAxes()
p3 <- FeaturePlot(object = Prog.data,
features = "Pseudotime",
pt.size = 2,
cols = rev(colorRampPalette(brewer.pal(n =10, name = "Spectral"))(100)),
reduction = "spring",
order = T) & NoAxes()
p1 + p2 + p3 Trajectories.progenitors <- Prog.data@meta.data %>%
select(Barcodes, nUMI, Spring_1, Spring_2, Pseudotime) %>%
mutate(Cycling.axis= Prog.data$Revelio.cc,
Phase = Prog.data$Revelio.phase,
Gmnc= Prog.data@assays$RNA@data["Gmnc",],
Ttr= Prog.data@assays$RNA@data["Ttr",],
Htr2c= Prog.data@assays$RNA@data["Htr2c",],
Top2a= Prog.data@assays$RNA@data["Top2a",])p1 <- ggplot(Trajectories.progenitors, aes(x= Pseudotime, y= Cycling.axis)) +
geom_point(aes(color= Phase), size=1.5) +
scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5]))
p2 <- Trajectories.progenitors %>% arrange(Gmnc) %>%
ggplot(aes(x= Pseudotime, y= Cycling.axis)) +
geom_point(aes(color=Gmnc), size=1.5) +
scale_color_gradientn(colours =c("grey90", brewer.pal(9,"YlGnBu")))
p3 <- Trajectories.progenitors %>% arrange(Ttr) %>%
ggplot(aes(x= Pseudotime, y= Cycling.axis)) +
geom_point(aes(color=Ttr), size=1.5) +
scale_color_gradientn(colours =c("grey90", brewer.pal(9,"YlGnBu")))
p4 <- Trajectories.progenitors %>% arrange(Htr2c) %>%
ggplot(aes(x= Pseudotime, y= Cycling.axis)) +
geom_point(aes(color=Htr2c), size=1.5) +
scale_color_gradientn(colours =c("grey90", brewer.pal(9,"YlGnBu")))
p1 + p2 + p3 + p4 + patchwork::plot_layout(ncol = 2)p1 <- ggplot(Trajectories.progenitors, aes(x= Cycling.axis, y= Gmnc)) +
geom_point(aes(color= Phase), size=0.5) +
scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) +
geom_smooth(method="loess", n= 50) +
ylim(0,NA)
p2 <- ggplot(Trajectories.progenitors, aes(x= Cycling.axis, y= Top2a)) +
geom_point(aes(color= Phase), size=0.5) +
scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) +
geom_smooth(method="loess", n= 50) +
ylim(0,NA)
p1 + p2rm(list = ls()[!ls() %in% c("Trajectories.progenitors", "ChP.data")])
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 5833769 311.6 10748997 574.1 10748997 574.1
## Vcells 59947732 457.4 391318918 2985.6 621902469 4744.8
Import progenitors cycling coordinates in the full dataset
ChP.data$Cycling.axis <- sapply(ChP.data$Barcodes,
FUN = function(x) {
if (x %in% Trajectories.progenitors$Barcodes) {
x = Trajectories.progenitors[x, "Cycling.axis"]
} else {
x = NA
}
})#date
format(Sys.time(), "%d %B, %Y, %H,%M")## [1] "09 December, 2021, 17,28"
#Packages used
sessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS/LAPACK: /shared/ifbstor1/software/miniconda/envs/r-4.1.0/lib/libopenblasp-r0.3.15.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] wesanderson_0.3.6 cowplot_1.1.1 ggExtra_0.9
## [4] ggplot2_3.3.5 RColorBrewer_1.1-2 dplyr_1.0.7
## [7] Matrix_1.3-4 orthologsBioMART_0.1.0 data.table_1.14.2
## [10] biomaRt_2.48.3 princurve_2.1.6 Revelio_0.1.0
## [13] SeuratObject_4.0.3 Seurat_4.0.5
##
## loaded via a namespace (and not attached):
## [1] BiocFileCache_2.0.0 plyr_1.8.6 igraph_1.2.9
## [4] lazyeval_0.2.2 splines_4.1.0 listenv_0.8.0
## [7] scattermore_0.7 GenomeInfoDb_1.30.0 digest_0.6.29
## [10] htmltools_0.5.2 fansi_0.5.0 magrittr_2.0.1
## [13] memoise_2.0.1 tensor_1.5 cluster_2.1.2
## [16] ROCR_1.0-11 globals_0.14.0 Biostrings_2.62.0
## [19] matrixStats_0.61.0 spatstat.sparse_2.0-0 prettyunits_1.1.1
## [22] colorspace_2.0-2 rappdirs_0.3.3 blob_1.2.2
## [25] ggrepel_0.9.1 xfun_0.28 crayon_1.4.2
## [28] RCurl_1.98-1.5 jsonlite_1.7.2 spatstat.data_2.1-0
## [31] survival_3.2-11 zoo_1.8-9 glue_1.5.1
## [34] polyclip_1.10-0 gtable_0.3.0 zlibbioc_1.40.0
## [37] XVector_0.34.0 leiden_0.3.9 future.apply_1.8.1
## [40] BiocGenerics_0.40.0 abind_1.4-5 scales_1.1.1
## [43] DBI_1.1.1 miniUI_0.1.1.1 Rcpp_1.0.7
## [46] viridisLite_0.4.0 xtable_1.8-4 progress_1.2.2
## [49] reticulate_1.22 spatstat.core_2.3-1 bit_4.0.4
## [52] stats4_4.1.0 htmlwidgets_1.5.4 httr_1.4.2
## [55] ellipsis_0.3.2 ica_1.0-2 farver_2.1.0
## [58] pkgconfig_2.0.3 XML_3.99-0.8 dbplyr_2.1.1
## [61] sass_0.4.0 uwot_0.1.10 deldir_1.0-6
## [64] utf8_1.2.2 labeling_0.4.2 tidyselect_1.1.1
## [67] rlang_0.4.12 reshape2_1.4.4 later_1.3.0
## [70] AnnotationDbi_1.54.1 munsell_0.5.0 tools_4.1.0
## [73] cachem_1.0.6 cli_3.1.0 generics_0.1.1
## [76] RSQLite_2.2.8 ggridges_0.5.3 evaluate_0.14
## [79] stringr_1.4.0 fastmap_1.1.0 yaml_2.2.1
## [82] goftest_1.2-3 knitr_1.36 bit64_4.0.5
## [85] fitdistrplus_1.1-6 purrr_0.3.4 RANN_2.6.1
## [88] KEGGREST_1.32.0 pbapply_1.5-0 future_1.23.0
## [91] nlme_3.1-152 mime_0.12 xml2_1.3.3
## [94] rstudioapi_0.13 compiler_4.1.0 filelock_1.0.2
## [97] curl_4.3.2 plotly_4.10.0 png_0.1-7
## [100] spatstat.utils_2.2-0 tibble_3.1.6 bslib_0.3.1
## [103] stringi_1.7.6 highr_0.9 lattice_0.20-44
## [106] vctrs_0.3.8 pillar_1.6.4 lifecycle_1.0.1
## [109] spatstat.geom_2.3-0 lmtest_0.9-39 jquerylib_0.1.4
## [112] RcppAnnoy_0.0.19 bitops_1.0-7 irlba_2.3.3
## [115] httpuv_1.6.3 patchwork_1.1.1 R6_2.5.1
## [118] promises_1.2.0.1 KernSmooth_2.23-20 gridExtra_2.3
## [121] IRanges_2.28.0 parallelly_1.29.0 codetools_0.2-18
## [124] MASS_7.3-54 assertthat_0.2.1 withr_2.4.3
## [127] sctransform_0.3.2 S4Vectors_0.32.0 GenomeInfoDbData_1.2.6
## [130] mgcv_1.8-36 parallel_4.1.0 hms_1.1.1
## [133] grid_4.1.0 rpart_4.1-15 tidyr_1.1.4
## [136] rmarkdown_2.11 Rtsne_0.15 Biobase_2.52.0
## [139] shiny_1.7.1
Institute of Psychiatry and Neuroscience of Paris, INSERM U1266, 75014, Paris, France, matthieu.moreau@inserm.fr↩︎